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Anisotropic classical Heisenberg models with all-to-all spin coupling display a Topological Non- 
connectivity Threshold (TNT) for any number iV of spins. Below this threshold, the energy surface 
is disconnected in two components with positive and negative total magnetizations respectively, so 
that magnetization cannot reverse its sign and ergodicity is broken, even at finite N. Here, we 
solve the model in the microcanonical ensemble, using a recently developed method based on large 
deviation techniques, and show that a phase transition is present at an energy higher than the 
TNT energy. In the energy range between the TNT energy and the phase transition, magnetization 
changes sign stochastically and its behavior can be fully characterized by an average magnetization 
reversal time. The time scale for magnetic reversal can be computed analytically, using statistical 
mechanics. Numerical simulations confirm this calculation and further show that the magnetic 
reversal time diverges with a power law at the TNT threshold, with a size dependent exponent. This 
exponent can be computed in the thermodynamic limit (N — > oo), by the knowledge of entropy as 
a function of magnetization, and turns out to be in reasonable agreement with finite N numerical 
simulations. We finally generalize our results to other models: Heisenberg chains with distance 
dependent coupling, small 3D clusters with nearest neighbor interactions, metastable states. We 
conjecture that the power-law divergence of the magnetic reversal time scale might be a universal 
signature of the presence of a TNT. 

PACS numbers: 05.45Pq, 05.45Mt, 03.67,Lx 



I. INTRODUCTION 

Statistical mechanics deals with systems containing a 
very large number (10 23 ) of interacting particles. Nowa- 
days, as the experimental investigation of few-atom sys- 
tems is becoming possible, the analysis of small systems 
raises fundamental questions and the problem of a 
statistical description of few-body systems with strong 
nonlinear interaction is a subject of current research 0. 
Unfortunately we are still far from understanding what 
are the conditions for a few-body system to reach, if any, 
an equilibrium, and how to describe it in the same way 
as statistical mechanics provides a powerful description 
of large systems. For instance, even the existence of tem- 
perature at the nanoscale has been recently questioned 
in Ref. 0. 

It has been recently shown that, for a Heisenberg 
model with all-to-all coupling, there exists a specific en- 
ergy threshold below which total magnetization cannot 
change its sign, even when the number of spins is fi- 
nite J4(. This ergodicity breaking phenomenon has been 
related to the existence of a Topological Non-connectivity 
Threshold (TNT) of the energy surface. This type of er- 
godicity breaking at finite N is different from the TV — > oo 
ergodicity breaking in a standard Ising model below the 
critical temperature. The existence of this threshold is 
not restricted to the infinite range coupling case. It is also 
present when the interaction among Heisenberg spins de- 
cays as R~ a , where R is the distance between two spins. 
It has been indeed proved [f| that, for a d— dimensional 



system, the ratio of the disconnected portion of the en- 
ergy range with respect to the total energy range tends 
to zero in the thermodynamic limit for a > d (short 
range interactions) while it remains finite for a < d (long 
range interactions). On the other hand, although the 
mean-field (all-to-all) type of spin coupling might ap- 
pear unphysical, magnetic systems can be realized, using 
modern experimental techniques [||, which are well de- 
scribed by Heisenberg-like Hamiltonians with an infinite 
range term. Moreover, when the range of the interaction 
is of the same order of the size of the system, all-to- 
all coupling may be a meaningful first order approxima- 
tion 0, Q ■ This could be the case for small systems used 
in current nano-technology, which requires to deal with 
systems made of a few dozens of particles. Otherwise, all- 
to-all coupling is relevant for macroscopic systems with 
long range interactions, like gravitational and unscreened 
Coulomb systems Q. 

We address in this paper the issue of providing a theo- 
retical framework to calculate the magnetization reversal 
time for the mean-field anisotropic Heisenberg model in 
a magnetic field, in which a finite number N of spins in- 
teract with all-to-all couplings. We have already stressed 
that below the TNT the total magnetization cannot re- 
verse its sign, thus magnetization does not relax to its 
equilibrium value. In this paper we will answer the fol- 
lowing questions: i) Above the TNT, does the magneti- 
zation reverse its sign? If so, on which time scale? ii) 
What is the relevance of the TNT in a system with a 
standard magnetic phase transition? Thus, the aim is 
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to explain the main physical effects associated with the 
TNT, and how it affects the phase transition appearing 
in this model at a higher energy. The latter is studied in 
the microcanonical ensemble, applying a recently devel- 
oped solution method of mean-field Hamiltonians based 
on large-deviation theory Q ■ We study in detail, numer- 
ically and analytically, the time scale for magnetization 
reversal. At the TNT, the reversal time diverges as a 
power law, with a characteristic exponent proportional 
to the number of spins N. Based on analytical calcula- 
tions, we expect this property to be universal. Finally, 
we extend some of the results obtained for all-to-all cou- 
pling to other models: chains with distance dependent 
couplings and small clusters with nearest neighbor inter- 
actions. We show the existence of the TNT also in these 
cases, and we present strong evidence for the power law 
divergence of the reversal time. 



II. THE MODEL 

The Hamiltonian of the model is 

N T N 



h = bY / S! + ^ X»< ' v : - ^sp, (i) 



,=1 



t=l j^i 



where Si = (Sf , Sf , Sf) is the spin vector with contin- 
uous components, N is the number of spins, B is the 
rescalcd external magnetic field strength and J the all- 
to-all coupling strength (the summation is extended over 
all pairs). Let us also define 
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as the components of the total magnetization of the sys- 
tem. Due to the anisotropy of the coupling, the system 
has an easy-axis of the magnetization along the y di- 
rection (the easy-axis of the magnetization is defined by 
the direction of the magnetization in the minimal energy 
configuration of the system). The equations of motion 
are derived in a standard way from Hamiltonian Q , and 
we obtain: 
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-BSf-JSfj: {j) S] 



BSf — JS'i 



(2) 



The total energy E = H and the spin moduli |5,| 2 = 1 
are constants of the motion. Dynamics, already stud- 
ied in a similar model 0, Q , is characterized by chaotic 
motion (positive maximal Lyapunov exponent) for not 
too small energy values and spin coupling constants. For 
J = the model is exactly intcgrable, while for generic 
J and B there is a mixed phase space with prevalently 
chaotic motion for \E\ < JN. 



III. THE TWO THRESHOLDS 

We will now show the existence of two distinct thresh- 
olds in this model: first we derive analytically the Topo- 
logical Non-connectivity Threshold (TNT), then we will 
present the microcanonical analysis and the analytical 
evaluation of the statistical threshold, at which a second 
order phase transition occurs in the N — > oo limit. 



A. The Topological Non— connectivity Threshold 

The phase space of the system is topologically discon- 
nected below a given energy density tdis, which can be 
obtained as in Ref. 0, . From symmetry considera- 
tions, both positive and negative regions of m y exist on 
the same energy surface. Indeed the Hamiltonian is in- 
variant under a rotation of 7r around the z axis for which 
Sf — > —Sf and Sf — > —Sf. Switching dynamically from 
a negative m y value to a positive one requires, for con- 
tinuity, to pass through m y = 0. Hence, for all energy 
values above 

e dls = Min[H/N \ m y = 0] 

magnetization reversal is possible, while below this value 
magnetization cannot change sign. 

Hamiltonian Q can be written as follows: 

H = BNm z + {ml - m 2 ) + I ^{Sff {Sf?. 

i 

(3) 

The Topological Non-connectivity Threshold (TNT) is 
defined as the minimum of the Hamiltonian under the 
N + 1 constraints: 

a) (Sf) 2 + {S?) 2 + (S?) 2 = l (4) 

b) m y = 0. (5) 

Instead of solving the constrained problem, we simplify 
it by calculating the absolute minimum of 

F = BNm z - J -Y}{Sff-{Sm 

If the minimal solution satisfies both m x = and m y = 0, 
the problem is equivalent to the original one. Conditions 
Q are taken into account setting: 

Sf = cos 9i , Sf = sin 9i cos fa , Sf = sin Bi sin fa . 

Taking the derivatives of F we obtain: 



dF . 2 

— — = J sin 0i cos d>i sin m = 
dF 



(6) 



^— = sm. 8 i[B + J cos 9i cos 2 4>i] = . (7) 

If B > J, Eq. has the solution, sin^ = 0, that also 
satisfies Eq. ©. It corresponds to all spins lying along 
the z-axis and 



-B 



(8) 
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If B < J, then from Eq. J7J) wo have two possible solu- 
tions for each i: 

1) Bi = tt: 

2) sin</>i = and cos 6*; = —B/J 

Let us define < n z < N as the number of spins satisfy- 
ing condition 1) above. Then 



F(n z ) = ^(B-J) 2 



N 



El 

2,1 



J 



so that the minimum is reached for n z = or cos 9i = 
—B/J and sin</>j = for all i. This in turn implies m y = 
and, for N even, m x = (choosing for instance 4>i = 
tt/2 for i = 1, N/2 and = -tt/2 for i = AT/2 + 1, TV). 
Then we have (for N even): 



BP_ 
2J 



Summarizing, we get [Tlj . 



-B for J < B 



-(§7 + D for 3 > B - 



(9) 



(10) 



The existence of e^is docs not represent a sufficient 
condition in order to demagnetize a sample for e > tdis- 
As it will be shown in Sec. II V HI regular structures in- 
deed appear in some cases, preventing most trajectories 
to cross the m y = plane. 



B. The Statistical Threshold: phase transition 

We now determine the statistical phase-transition en- 
ergy of the model in the microcanonical ensemble. To 
keep the calculations easy, wc will first neglect the term 



y\2 



take into account this term, 
calculations, we will also set 



{Sf) 2 in ©. We will show later how to 
In order to facilitate the 



e/B 
JN 
~B 



(11) 



Thus we can consider the following Mean-Field Hamil- 
tonian: 



H 



MF 



N 



(12) 



Note that this mean field limit is formally identical to 
phenomenological single spin Hamiltonians used to model 
micromagnetic systems [l2|. 

Using this simplified Hamiltonian, we can calculate the 
entropy, counting the number of microscopic configura- 
tions associated with given values of m„ m y and m z , 



independently of the energy of the system. This can be 
done using Cramer theorem, a basic tool of Large Devi- 
ation Theory 0|. Each single spin is characterized by 
two angles 9 and (j>, such that S z = cos 6, S x = sin 9 cos </>, 
sin 6. We calculate the function 



S y = sin 



— / sm9 d9 / dd> 
^ h Jo 

A a; sin 8 cos 4>-\-\ y sin 6 sin <f> 



(13) 



We then get the entropy s(m x ,m y ,m z ) through a 
Legendre-Fenchel transform of ln\E': 



i(m x ,m y7 m z ) = - sup [X x m x 



XyUly 



Am, 



ln*(A x ,A y ,A)] 



(14) 



This calculation gives us an approximate expression for 
the probability P(m x , m y , e), which describes the system 
for each energy: 



P(m x , m y , e) oc exp [Ns(m x , rn y , \ 



K))) 

(15) 

Integrating over m x , one gets the marginal probabil- 
ity distribution P(m y ,e). We define the paramagnetic 
(resp. ferromagnetic) phase by a probability distribution 
P{m y ,e) which is single peaked around m y = (resp. 
double peaked). To locate the statistical phase transi- 
tion energy e sta t, we assume that the transition is second 
order; it is then sufficient to study the entropy around 
iriy = 0. We will also set m x = 0, since it is easy to see 
that a non-zero m x would only decrease the entropy for 
negative energy states; thus these states with non zero 
m x have little influence. Physically, the picture is the 
following: the negative energy has to be absorbed by ei- 
ther a non-zero m z , or a non-zero m y or both. For small 
negative energies, it is entropically favorable to decrease 
a bit ni z , since it has a linear effect on the energy. For 
negative enough energies however, it costs much entropy 



to decrease 



further, so that a non-zero : 



is favored, 



this is the phase transition. As a small m y results in a 
small A y , we develop and In 'J up to second order in 
A„ : 



*(A„A) 
ln*(A a ,A) 



sinh A \ y A cosh A 



sinh A 



A 

. / sinh A 



A 3 

X y A cosh A 



(16) 



sinh A 



A 2 sinh A 



(17) 



The maximization over A and X y yields the equations: 



m z - A0(A) + ^0'(A) 

m y = Xy(f){\) , 



where 



A cosh A — sinh A 
A 2 sinh A 



(18) 
(19) 

(20) 
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From Eq. JTSJ, wc write A = Ao + 02 A^, where Ao is 
defined implicitly by 



m z = A o 0(A o ) 



(21) 



and 02 is a coefficient. We then compute the entropy up 
to second order in m y : 

s(m y ,e) = - 



m v > , \ i { sinhAn , 

-A (m z )+ln — ^ ; (22) 



2<MA) 



V A 



note that the terms with canceled. Using energy con- 
servation m z = e+Irriy/2, we obtain the entropy s(m y ; e) 
as a function of m y alone, e being now a parameter. The 
equation for Aq is: 



e+ 2 m ? = A o<MA ) 



(23) 



We then write Ao = \i + m 2 ^ 2 , with e = /x(/)(/i), and 
substitute into Eq. i|22|) . to get s(m y ;e) up to second 
order in m„: 



s(m y , 1 



-/ie + In 



/ sinh 



to,, 



1 

2e 



(24) 



The vanishing of the second derivative in m y = yields 
the critical energy: e s t a t = — 1/-^ which can be expressed 
in the old variables, see Eq. I|ll() : 
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FIG. 1: Comparison of the distributions P(m y ) obtained ana- 
lytically (solid lines) and numerically (open circles). All plots 
are for N — 6 spins, a coupling J — 1/3 and a field B — 1. 
From top to bottom, the energy per spin is e = —0.7 (fer- 
romagnetic phase), e — —0.5 (close to the phase transition), 
and e = —0.3 (paramagnetic phase). In panel a) the dashed 
line shows the improved analytical calculations in which the 
non-mean field terms are taken into account. 



B 2 
JN' 



estat = "Tat- ( 2 ^) 



At this threshold entropy has a maximum in m y = 0, 
with vanishing second derivative. In the thermodynamic 
limit the second derivative of the entropy as a function of 
e becomes discontinuous in e stat , indicating that a true 
second order phase transition occurs at e s t a t-, for N — > 00. 
This analytically calculated value of e stat is in reasonable 
agreement with numerical results obtained using the full 
Hamiltonian 0J. 

The corresponding probability distribution P(m y ,e) 
obtained from the mean field Hamiltonian (|12f) . should 
be compared with that obtained (numerically) from 
the full Hamiltonian QJ, for instance by sampling of 
the phase space. Results are shown in Fig. As one 
can see the agreement is quantitatively good in the 
paramagnetic phase, but only qualitatively correct in 
the ferromagnetic phase; here the double peaked shape 
is correct, but the details are significantly off. 

The inaccuracy of the calculation may come from both 
the small value of N, and from the term J /2}] ! (S'f) 2 — 
(Sf) 2 , which has been neglected till now. It can be in- 
cluded in the statistical analysis as follows. The Hamil- 
tonian depends now on another global quantity, A = 
((Sf) 2 — (Sf) 2 }. It is possible to include it in the large 
deviation calculation; Eq. (|13fl is modified into: 



V(\ X , Xy,X, /J.) 



1 

4~ 



sin 9 d6 

Jo 



gAa; sin 9 cos 0+A y sin 9 sin ^sin 0(cos 0— sin tf>) (26) 



One proceeds by writing a probability distribu- 
tion P(m x , m y , m z , A) cx exp(7Vs(m 2: , m y , m z , A)), 
taking into account the energy conservation 
m z = e - I (ml - m 2 )/2 + JA/2 (for B = 1), and 
integrating over m x and A to obtain P(m yi e). This 
last step has to be carried out numerically, and no 
simple expression as (|25|l is available any more. A 
comparison with a numerical investigation of the phase 
space shows that the additional term has a significant 
contribution; the P(m y ) we obtained in the ferromag- 
netic phase improves on the mean field calculation, 
see Fig. ^ We conclude that the remaining discrep- 
ancies come from the small value of N (N = 6 on Fig.^l. 

Let us finally remark that in this statistical framework, 
the TNT energy, e^s , can be recovered as the energy such 
that s(0,e dls ) = — oo. From Eqs. (|21I24() with m y = 0, 
it is easy to get Cdis = — 1; this implies, using Eq. ifTTjl. 
(-dis = — B, which is the same as in Eq. l|lUf) in the limit 
N -> oo for J < B. 



IV. TIME SCALE FOR MAGNETIC REVERSAL 

In the following, we will study the dynamics of the full 
Hamiltonian iJTjl. which, at variance with (U2J, is non- 
integrable and can display chaotic motion. Let us first 
notice that in the large N limit the minimal energy can 
be easily estimated (see Appendix as 
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FIG. 2: In this figure, all data refer to the N = 6, B = 1, 
J = 3 case, a) Magnetization m y vs time. Reversal times 
(times between neighbor zeros of m y have been indicated as 
n, . . . , Tfc. b) Solid line : probability distribution of reversal 
times versus the normalized time r/(r), for 1 trajectory and 
10 4 different crossings. Dashed line : probability distribution 



of relaxation times P — 



p + (0)-i/2 vs 2i/( T )- As initial condi- 
tion we choose P+ (0) = 1 and 10 4 different initial conditions. 



2JN 2 iul 



/ > 2- 

— N 



(27) 



-B 



for 



/ < - 

AT" 



From Eqs. l)27l2^ll0|) we have that il J > B/N then 
fstat > £<iis > £?rim- In what follows we will restrict our 
consideration to the region of parameters for which these 
three thresholds are different. 

The two thresholds, edis and e s tat, define three energy 
regions which show different dynamical and statistical 
properties: 

1) For e < e,us, the probability distribution of m y P(m y ) 
has two separate peaks, with P(m y = 0) = 0, so that m y 
cannot change sign in time. 

2) For e s tat < e < 0, m y quickly changes sign in time and 
P(m y ) is peaked at m y = 0. 

3) For edis < e < e s tat, the probability distribution is 
doubly peaked around the most probable values of the 
magnetization. These two peaks are not separated and 
P = P{ m y = 0) 7^ 0. What actually happens dynami- 
cally depends on the relative strength of the coupling J 
with respect to B. More specifically we can characterize 
two different behaviors, chaotic and quasi-intcgrablc. 



A. Chaotic Regime 

1. Time scale for magnetic reversal and relaxation 

For J big enough (fully chaotic regime) the behav- 
ior of m y (t) resembles a random telegraph noise |l5j| . 
Fig-Et) : magnetization switches stochastically between 
its two most probable values, reversing its sign at random 
times. If we sample the magnetization reversal times, 
Tfc, defined as the time interval between two crossings of 
m y — 0, we find that they follow a Poissonian distribu- 
tion with average (r) . Such distribution of the reversal 
times is a consequence of strong chaos: the system looses 
its memory due to sensitivity to initial conditions and 
the reversal probability per unit time. A = l/(r), be- 
comes time independent. 

Since the magnetization reverses its sign randomly, any 
initial macroscopic sample with m y 0, will relax to an 
equilibrium distribution with a vanishing average mag- 
netization. In order to characterize quantitatively the 
relaxation process, we introduce the probability to have 
a positive magnetization, P+(i), at time t. This is mea- 
sured by considering an ensemble of n initial conditions 
and counting, for each time t, the number of trajecto- 
ries n + (t) for which m y > 0. At equilibrium we have 
P + = 1/2, in agreement with standard statistical me- 
chanics considerations. Below edis, P+ cannot change in 
time because the sign of m y remains the same for all 
trajectories. Above edis, P+{t) can change in time. Nu- 
merical results show that P+(t) decays exponentially to 
the equilibrium value 1/2 and that the time scale for 
reaching the equilibrium value is independent of the ini- 
tial probability distribution, P + (0), see Fig.[2)) (dashed 
line) . 

A simple statistical model can explain the qualitative 
features of this magnetic relaxation process. Let us start 
with an ensemble of n initial conditions, of which n+ 
with a positive magnetization and n_ with a negative 
magnetization, such that n = n + + n_ . Assuming that 
m y can take only two values, + and — , we can write 
a pair of differential equation for the populations with 
positive and negative magnetizations: 

h + = — Xn + + An_ 
n_ = — An_ + An+. 

Where A is the reversal probability per unit time defined 
above. Defining P + = n+/n, we can solve these equa- 
tions, obtaining 



P+{t) 



1 



(P + (0)--)e 



-2Xt 



(28) 



P+(t) reaches the equilibrium value with a typical re- 
laxation time t = 1/(2A). This simple model predicts 
a magnetic relaxation time, r, proportional to the av- 
erage magnetization reversal time, which checks pretty 
well with numerics (compare solid and dashed lines in 
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FIG. 3: Panel a) shows the scaling of the magnetic relaxation 
time t for different values of N, for the case B = 0, J = 3. 
The relaxation time has been computed starting from an en- 
semble of 100 different initial conditions. The dashed line is 
the linear fit: -^log{r) = 0.075 — 0.82log(e — e<jj s ). In panel 
b) we show the times resulting from statistical mechanics, see 
Eq. 131H . vs the dynamically computed relaxation times, for 
different N values. Pmax/Po, the main ingredient of formula 
131(1 is numerically determined for each energy density choos- 
ing 10 6 points in a small energy interval. Symbols for the sta- 
tistical times are as follows: N — 6 (circle), N = 12 (cross), 
N = 24 (plus), N = 48 (diamond). Dashed line is y = x. 



Fig. EJd). Therefore, hereafter we will use indifferently 
the two concepts. 

Analyzing the magnetic relaxation times for all ener- 
gies in the range {eais, tstat), we find that they grow expo- 
nentially with the number of spins for sufficiently large A, 
as expected for mean-field models. More remarkable is 
the power law divergence of relaxation time at the non- 
connectivity threshold. Numerical data are consistent 
with the following scaling law, 



(29) 



for which a theoretical justification will be given be- 
low. Eq. I|29|) is valid above the non-connectivity thresh- 
old and not too close to the statistical threshold e s tat- 
The comparison of this formula with numerical results is 
shown in Fig. |3J 

To explain and substantiate these numerical findings, 
we now turn to an analytical estimate of the relaxation 
times, based on statistical mechanics. In Refs. [Tti Il7j. 
on the basis of fluctuation theory [l8L Il9j , it has been 
argued that metastable states relax to the most proba- 
ble state on times proportional to exp(AAs) where N 
is the number of degrees of freedom and As is the spe- 
cific entropic barrier. In our case exp(AAs) is nothing 
but Pmax/Po, where P max is the value of P{m y ) for the 
most probable value of m y , and Pq = P(m y = 0). Thus, 



the exponential divergence as a function of A shown in 
Fig. Eh) is consistent with Refs. [l6l fl7l | . These papers, 
however, did not study the behavior of t at fixed A in 
the neighborhood of the non-connectivity threshold. We 
perform this calculation in Appendix iBl obtaining 



V( e_ 6 dis) 



(30) 



with a — 1 gcnerically, but a = 3/4 for B = J = 1. 
This result is qualitatively correct (power law divergence, 
exponent proportional to N), and quantitatively reason- 
able. Indeed, numerical simulations give a ~ 0.82 (in- 
stead of a = 1) for J 3> B and B = (see Fig. and 
a ~ 0.55 (instead of a = 3/4) for B = J = 1. We expect 
these qualitative features to be valid beyond the all-to-all 
coupling studied here, as it will be shown in Sec. [VJ 

The calculations to evaluate P ma x/Po rely on several 
approximations, the most doubtful being the large N as- 

Hence, 



sumption (as seen also in Section IlIIBIl . Hence, despite 
the discrepancies in the exponents found above, the pro- 
portionality between r and P ma xl Po may still be valid, 
also for small N. To test this proportionality, we have 
calculated numerically the value of P ma x I Po i an d we have 
found this value to be proportional to the relaxation time 
in any case. In particular, for the case B = 0, we have 
found a very good fit setting 



2 P 



J Po 



(31) 



The Pmax/Po factor in this formula represents the prob- 
ability to cross the entropic barrier, and the 1/J factor 
can be heuristically associated with the typical time scale 
of the system (for B = the Hamiltonian is proportional 
to J). A deeper theoretical justification of this formula 
should be obtained in view of its success in describing the 
numerical results for different TV values (see Fig.©)). 



2. Chaotic Driven Phase Transition 

Let us now answer the following question: if the 
measured values of the magnetization are given by the 
time average of the magnetization, for which energies 
the system will be found magnetized and for which 
unmagnctized? 

From Eq. l|29|l and from the proportionality of the re- 
laxation times and the reversal times, it is clear that the 
infinite time average of the magnetization will be zero 
above the TNT and different from zero below, due to the 
divergence of the reversal time. Nevertheless, the con- 
clusion is different for a finite observational time t i> s . In 
Fig. 0] we show the time- averaged magnetization 

1 r^' 
{m y )obs = / dt m y (t) 

Tabs Jo 

vs the specific energy e for N = 5 (Fig. HJi) an d A = 50 
(Fig. 0})) spins during a fixed observational time. While 
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FIG. 4: Time average of m y over the observational time Tobs 
vs e for different number of particles (a) TV = 5, (b) TV = 50, 
with fixed J — B — 1. Each single point has been obtained 
taking the time average over the time intervals r b 3 = 10° 
(a) and r b s = 10 4 (b). Dashed curves indicate the equilib- 
rium value of m y , obtained from statistical mechanics. Ver- 
tical lines represent the non-connectivity and the statistical 
threshold respectively. The arrow in panel (b) indicates the 
energy value e bs of the chaotic driven phase transition due 
to the finite observational time. 



in (a) (m y ) b a is zero just above e d i S , in (b) it vanishes 
at a value e b s located between edis and e s t a t- Indeed, if 
Tobs ^ t, the magnetization has time to flip between the 
two opposite states and, as a consequence, (m y ) b s — 0. 
On the contrary, if r b s <C r the magnetization keeps its 
sign and cannot vanish during r bs . Defining an effective 
transition energy e b s from r b s = T(e bs), one gets, in- 
verting Eq. (|29|l . the value indicated by the vertical arrow 
in Fig. |^). This is, a posteriori, a further demonstration 
of the validity of Eq. for any N. 

From a theoretical point of view, it is interesting to 
note that, for any fixed N, if the fully chaotic regime 
persists down to e dis , e obs — ► £dis when T obs — ► oo. On 
the other hand, in agreement with statistical mechanics, 
for any finite T obs , e a bs — ► e s tat when N — > oo. This 
implies that the limits r b s — * oo and N — > oo do not 
commute. From the above considerations it follows 
that if T bs — > oo at finite N, the threshold which 
distinguishes between a magnetized energy region and 
an unmagnctized one is edis and not e sta t- We can 
thus consider edis as the critical threshold at which a 
"dynamical" phase transition takes place: we call this 
transition a chaotic driven phase transition. 

Let us finally note that, usually, for long-range inter- 
actions, the interaction strength is rescaled in order to 
keep energy extensive |2(J . In our case this can be done 
setting J = I/N. With this choice of J as N — > oo, at 
fixed I, J becomes much smaller than B, then a quasi- 
integrable regime sets in and Eq. I|29|) looses its validity 



(see Sec. HVBl 
hidden. 



The presence of the TNT is therefore 



B. Quasi— integrable Regime 

In this Section we will give numerical evidence of the 
quasi-integrable regime for J < B, in the energy region 
between e d i S and e s t a t- If the system dynamics is not in 
a fully chaotic regime, there are important consequences 
for reversal times. For instance reversal times strongly 
depend on initial conditions and Eq. i|29|) looses its va- 
lidity: 

In Fig.[S]we consider a system with different interaction 
strengths J in order to enter a quasi-integrable regime (a 
and b) and a chaotic one (c and d). The energy in the 
two cases has been chosen such that the entropic barrier 
is roughly the same, see Fig. [JJ b) and d). This means 
that, from a statistical point of view, both systems are 
characterized by roughly the same probability to jump 
over the barrier. Nevertheless, as one can see in pan- 
els a) and c) of Fig. [SJ the behavior of the probability 
^(*) = p+(o)-!/2 significantly differs in the two cases. 
Such a big difference in the statistical properties of mag- 
netic reversal times can be explained only by a drastic 
change in the dynamical properties of the system. In- 
deed, while Fig. UJi refers to a quasi-integrable regime, 
Fig. |5j; refers to a fully chaotic regime. This cannot be 
explained by the different J values, which, as we have 
shown in the previous Section, have only a linear effect 
on the reversal probability per unit time. 

To better understand the origin of this quasi-integrable 
regime, it is interesting to compare the dynamics ob- 
tained from the full Hamiltonian J3Jl with the dynamics 
obtained from the Mean-Field Hamiltonian (|I2|I . Tak- 
ing into account the conservation of the total angular 



momentum M 2 = m 2 



a change of variable 



maps l)12|l onto a 2 degrees of freedom Hamiltonian; the 
dynamics of the global magnetization is then obviously 
integrable. 

In Fig. we show the projection of some trajectories 
on the (m x ,m y ) plane. We considered the two different 
dynamical regimes described above. For definitcness, we 
vary J but we choose the specific energy in order to keep 
the same value of P ma x/Po ~ 20. Let us first discuss 
Fig. Dark lines represent orbits of the Mean-Field 
Hamiltonian Q12|). The orbits of the macroscopic vari- 
able fh = (m x ,m y ,m z ) cover tori, since the Mean-Field 
Hamiltonian l|12|l is exactly integrable. Nevertheless, tra- 
jectories display different features: while trajectory (I) 
crosses the line m y = 0, trajectory (2) remains confined 
in the negative (m y < 0) branch belonging to the same 
energy surface. Two trajectories of the full Hamiltonian 
© and the same initial conditions are then considered, 
labeled by (3) and (4). As one can see these orbits stay 
for a long time sufficiently close to the Mean-Field or- 
bits. Again, while (4) displays a typical "ferromagnetic" 
behavior, (3) is of "paramagnetic" nature. Both trajec- 
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FIG. 5: (Color online) Magnetic relaxation probability P(t) — 
p + (0)-i/2 anc ' probability distributions of m y in the quasi- 
integrable (J = 0.5, e = —0.8) (red panels a),b)), and chaotic 
(J — 3,e = —0.9) (black panels c),d)) regimes. For both cases, 
TV = 6 and B = 1. The specific energies for the two cases 
have been chosen in such a way that Pmax/Po is the same 
(see panels b) and d)). We compare P(t) for the two different 
dynamical regimes. In the fully chaotic regime, panel c), P(t) 
decays exponentially with time t, and the average relaxation 
time is of order 1. In the quasi-integrable regime, panel a), 
the decay of P(t) is much slower (observe the difference in the 
time axis scale). 
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FIG. 6: (Color online) Projections of the trajectories of 
m = (jn x ,m y ,m~) over the (m x ,m y ) plane for both the 
Mean-Field Hamiltonian 21 (heavy full lines : (1) blue and 
(2) black ) and the full Hamiltonian @ (chaotically scat- 
tered dots (3) yellow and (4) green). In panel a) (J = 0.2, 
e = —0.3), the system is in a quasi-integrable regime. The 
orbits of the full Hamiltonian remain close to the orbits of 
the Mean-Field Hamiltonian for all the integration time con- 
sidered (10 4 ). Note also that in the some of the orbits cross 
the m v = line, thus demagnetizing the system, and some do 
not. In panel 6) ( J = 3, e = —0.32) the system is in a highly 
chaotic regime. In this case the orbit of the full Hamiltonian 
(6-yellow) does not remain close to the one of the Mean-Field 
Hamiltonian (5-black) and covers most of the available phase- 
space. The integration time is, also in this case, 10 4 . 



tories (3) and (4) have a positive maximal Lyapunov and 
are therefore chaotic. Upon increasing J, and keeping 
the same value of P m ax/Po ~ 20, we enter in the regime 
described by the lower panel in Fig. |S| In this case, as 
above, we show the orbit (5) of the Mean-Field Hamilto- 
nian (actually a "ferromagnetic" one). The correspond- 
ing orbit of the full Hamiltonian, (6), is still character- 
ized by a positive Lyapunov exponent and covers both 
branches, m y > and m y < 0, thus inducing the demag- 
netization of the system. What is important to stress is 
that in this case the trajectories of the full Hamiltonian 
cover both the positive and the negative magnetization 
branch on the same energy surface. Having in mind the 
mechanism that produces the transition to global stochas- 
ticity in low-dimensional Hamiltonian systems plj . we 
can conjecture that invariant curves, confining the mo- 
tion, exist in the case of Fig. [SJt. The breakdown of 
these invariant curves signals the transition to a globally 
chaotic motion. Of course, characterizing such a break- 
down is an hard task, due to the high-dimensionality of 
the phase-space. 

The determination of parameter regions in which the 
system is quasi-integrable or fully chaotic is still an open 
question. We can only make a few qualitative consider- 
ations. Let us consider Hamiltonian ©; it contains the 
sum of two terms: a mean-field integrable term plus the 
term J/2J2i(Sf) 2 ~ ( S"f ) 2 , which is responsible for the 



chaoticity of the system. The minimal specific energy of 
this term is e c ^ aos ~ —J/2. We can thus suppose that 
for e < e c haos the quasi-integrable regime prevails, while 
for e > e c haos a fully chaotic regime sets in. Thus, in 
order to have a fully chaotic regime in the energy region 
between e dis and e sta t, it is necessary that e dis > e chaos . 
This is always the case if J > 2B since for these values of 
J, e dis ~ e chaos- On the contrary, for J < 2B we expect 
a quasi-integrable regime between edis and e c haos, which 
should persist in the thermodynamic limit. 

V. OTHER MODELS 

Till now, we have concentrated our analysis on a spin 
system with all-to-all anisotropic coupling. The results 
obtained concerning the TNT and the time scales for 
magnetic reversal can be extended to more general situ- 
ations. In this Section, we consider two possible gener- 
alizations, and discuss how our results can be extended 
to: i) distance dependent interactions and ii) metastable 
states. 

i) Distance dependent coupling. In Ref.|j| it has been 
considered a spin coupling which decays with the dis- 
tance as R~ a . It is possible to prove that in the N — > oo 
limit, for a < d, a finite portion r of the energy range cor- 
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FIG. 7: Magnetization relaxation time r vs. energy density 
for two different models, (open symbols) a chain of N = 6 
Heisenberg spins with R~ a interaction and a — 2. A best fit 
gives log(r) ~ -3.9 log(e - edis), with e dis ~ -0.71. Here, 
E m ,„ = —1.083. (full symbols) 3-D cube with an additional 
spin at the center and nearest neighbor interaction. A best 
fit gives log(r) ~ —5.5 log(e — (. dis ), with e dis ~ —1.33. Here, 

£min — 20/9. 



responds to a disconnected energy surface. For a > d, 
this portion vanishes in the A — > oo limit. For finite 
N, however, a well defined non connectivity threshold 
£dis > £mm exists in both the short and the long case 
when the anisotropy of the coupling induces an easy-axis 
of the magnetization. Numerical simulations support the 
conjecture that the behavior of the average magnetiza- 
tion reversal time is qualitatively similar to the a = 
case. A power law divergence of the average reversal 
time when e approaches edis is observed, see Fig. {7\ 

More realistic models of micromagnctic systems in- 
clude 3D clusters of spins interacting only with their 
neighbors. Again, for large N, the non connectivity 
threshold energy edis converges to the ground state en- 
ergy e m i n . However, for small clusters a significant por- 
tion of the energy range corresponds to a disconnected 
energy surface. We have performed numerical simula- 
tions on a cluster of 9 spins, arranged on a cube, with one 
spin in the middle. Each spin of the cube interacts with 
its 3 neighbors and with the middle spin. Fig. [7| shows 
that the divergence of the magnetic relaxation time close 
to edis is again compatible with a power law. 

ii) Metastable states. The existence of the TNT 
has important consequences for the decay time from 
metastable states. In order to discuss this feature for a 
simple example, let us consider Hamiltonian QJ, adding a 
term, B y ^ Sf , which contains a coupling to an external 
field directed along the easy-axis of the magnetization. 
In this case the non-connectivity threshold still exists 
(it has the same value as before) but the two peaks of 
P(m v ,e) below e s tat do not have the same height, see 
Fig. Thus, we can consider the time needed to reach 
the equilibrium value of the magnetization if we start 
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FIG. 8: Decay time from a metastable state. In panel a) 
P(m y ) is shown for e = —0.5. In panel b) we show the power 
law divergence of the decay time from a metastable state (full 
circles). To compute this decay time, n = 100 initial condi- 
tions have been considered for each energy. The best linear 
fit log(r) = 0.309 — 5.27log(e — edis) is also shown (dashed 
line). The decay times computed from Eq. (1321 are shown 
as crosses. Apart from a deviation for high specific energy, 
the statistically computed decay times are in good agreement 
with the numerically computed ones (full circles). Parameter 
values are: N = 6, J = 3, B = 0, B y — 1 



from a metastable state. Below edis metastable states 
becomes stable for any finite N. Above edis the decay 
time diverges at edis as a power law, see Fig. |S|d. This 
decay time can be estimated from the statistical proper- 
ties of the system. Indeed, employing the same simple 
model described in Sec. IjlV A 1(1 . we can evaluate the de- 
cay time scale. Denoting by P^ax? ^max an d Pmin the 
probabilities of the thermodynamic stable, metastable 
and unstable state, respectively (see Fig. |HK), and set- 
ting A± = P± ax /P min , and q = P+ ax /P~ ax we get the 
following estimate of the decay time: 



1 Prr 



1 + qPn 



(32) 



The good agreement of this estimate with the computed 
decay times is shown by the crosses in Fig. OJ). 



VI. CONCLUSIONS 

Anisotropic Heisenberg spin models with all-to-all cou- 
pling show a Topological Non-connectivity Threshold 
(TNT) energy Below this threshold the energy sur- 
face splits in two components, with opposite easy-axis 
magnetizations, and ergodicity is broken, even with a fi- 
nite number of spins. In the same model, a second or- 
der phase transition is present, at an energy higher than 
the TNT energy. We have fully characterized this phase 
transition in the microcanonical ensemble, using a newly 
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developed method, based on large deviation theory [8(. 
For energies in the range between the TNT and the phase 
transition, magnetization randomly flips if certain strong 
chaotic motion features are present: the statistics of mag- 
netization reversals is Poissonian. Based on the knowl- 
edge of the microcanonical entropy as a function of both 
energy and magnetization, we have derived a formula for 
the average magnetic reversal time, which is valid in the 
large N limit. This formula agrees well with numerical 
results. The formula also predicts a power law divergence 
of the mean reversal time at the TNT energy, which is 
also well verified in numerical experiments. The expo- 
nent of the power-law divergence is also in reasonable 
agreement, although finite N effects are quantitatively 
important. 

Finally, we have shown that all these features (presence 
of a TNT, power-law divergence of the reversal time, etc.) 
are not limited to systems with all-to-all coupling. The 
phenomenology is qualitatively the same for anisotropic 
Heisenberg spin models with distance dependent inter- 
actions and for small clusters of Heisenberg spins with 
nearest neighbor coupling. We also considered systems 
where metastable states are present. In this case, while 
below the TNT they are trapped, above it their decay 
time diverges as a power law at the TNT. Therefore, we 
conjecture that the power-law divergence of the magnetic 
reversal time may be a universal signature of the pres- 
ence of the TNT, which is a generic feature of systems 
with long-range interactions or small systems for which 
the range of the interaction is of the order of system size, 
if the anisotropy of the coupling is such to determine an 
easy-axis of the magnetization. 
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and verify that it satisfies m x — 0. Taking derivatives 
dNe 



d(j)i 

dNe 
~d0~ 



Im y sin 0i cos <f>i = 



sin 6i — Im y cos 9i sin ( 



(A2) 



one gets two kinds of solutions (both with m x = 0): 

1) 0i = 7r and 4>i = 0, tt 

2) (pi = ±7r/2 and tan#, = ±Im y . 

Let us define Nn\ the number of solutions of type 1) and 
Nn 2 the number of solutions of type 2) so that ni+n 2 = 
1. Since m z = — m — n 2 cos 9 and m y = ±n 2 sin# where 
9 is the solution of type 2), condition 2) is equivalent to 
cos 9 = 1/In 2 . Therefore, when In 2 < 1 the set defined 
from 2) is empty and only solutions in the class 1) can 
be obtained. It is also easy to find the expression for the 
energy in terms of 1/1 < n 2 < 1: 



1 

e = - 1 -27 



712 - 2 n 2- 



(A3) 



Minima must be sought among the extrema so that when 
n 2 = 1 then e m i n = —1/2/ — 1/2 and when n 2 — 1/1 then 
e?nm = — 1- In terms of /, one then has 



-1/2/ -1/2 for J>1 



-1 



for 



(A4) 



/ < 1. 



From Eq. 

CD- 



4j) we have H10|) . using transformations in 
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APPENDIX A: MINIMUM ENERGY 

In this section we find the minimum of the Mean-Field 
Hamiltonian i|12fl : 



e = m z + -{m 2 x -m y ). 



It is sufficient to find the absolute minimum of 



(Al) 



m z - (I/2)m 2 y , 



In this section, we study the divergence of the reversal 
time for e — > ej is , at fixed N. Let us assume that it is 
given by 

r ~ maxP(m„, e)/P(m y = 0, e) , 
we show that 

(Bl) 



with a a constant independent of JV; we find a = 1 or 
a = 3/4, depending on the parameters of the Hamilto- 
nian. 

First, we note that although max m P(m y , e) increases 
exponentially with N at fixed e, it does not change much 



at fixed N when e 



the behavior of r is dominated 



by the value of P(m y = 0, e) cx e -^(™«=o,<0. The prob- 
lem is then reduced to the computation of s(m y = 0, e). 
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Before turning to the actual calculation of s(m y = 
0, e — ► edis), we consider the following problem, which 
will be useful later. Let us consider the random vari- 
able x in [—1,1] with distribution p(x); we call m = 
y~li—i Xj/N , and ask the question: what is the behavior 
of P(m = 1 — 8), 5 << 1, for fixed N reasonably large? 
Using Cramer's theorem, we write 



tt(A) 



p(x) exp(Xx)dx 



and 



s(m) = — sup (Am — In ^(A)) 

A 



(B2) 



(B3) 



m — > 1 implies A — > oo, so the behavior of p close to 
x = 1 dominates i|B2|) . We write ~ a(l — a;) 7 , close 
to x = 1, with 7 > —1. Then for A — » oo, 



(A) - a- 



Pe- U du 



Then the maximizing A in (|B3() is given by 

m = 1 — (7 + 1)/A. Substituting into (IB3I) . 

we get s(m = 1 — 8) ~ (7 + 1) ln<5, and finally, 
P(m = l-,5)-<5 Ar ( 7+1 ). 

We now apply this result to the easiest case, the sim- 
plified Hamiltonian H = N(m z + I(m 2 x - m 2 y )/2). The 
threshold is ejis = — 1 (we consider the case / > 1); we 
write e = —1 + 8, with 8 « 1. We want to compute the 
entropy s(m x ,m y = 0,m z = —1 + 8 — Ira 2 ,). Noticing 
that a small <5 implies a small m x , we simplify the calcu- 
lation to s(m x = 0, m y = 0, m z = —1 + 8). We now use 
once again Cramer's theorem. 

■1 pit /°2tt 

ty(\ x , A y , X z ) = — / sin9d9 / c?0cxp(A 2; sin cos0) 
47r Jo Jo 

exp(A y sin6* sin0 + A z cos6*) . (B4) 
Then s is given by: 



s(m x ,m y ,m z ) = - sup [m x X x + m y \ y + m z \ z 



- ^^(Aa, X y , A 2 )] 



(B5) 



For m x = m y = 0, the maximizing X x and A y are found 
to vanish. Thus, the problem reduces to calculating 
s(m z = 1 — 8) (using also the symmetry m z — ► —m z ). 
Recalling that m z = (cos 9), with 9 the latitude of 
a point taken randomly on the sphere with uniform 
probability. This is equivalent to saying that m z = (q), 
with q a random variable uniformly distributed between 
— 1 and 1. Using the general result derived above with 
7 = 0, we find that r ~ (e — e<iis) aN , with a = 1. 

We turn now to the complete Hamiltonian, with 
B = 0, H = N{I{m 2 x - ml)/2 - J A/2). The threshold 
is now Cdis = —J/2. We set e = —J/2 (1 — 8). Noticing 



that again, a small 8 implies a small m x , we compute 
s(m x = 0,m y = 0, A = 1 — 8), in the limit of small 8. 
A is defined as (sin 2 9 cos 2(f>) , for 9 and <fi coordinates 
of points taken randomly on the sphere with uniform 
probability. Again, this is equivalent to saying that 
A = (q), with q now having a non uniform distribution 
p(q) in [—1,1]. However, p(q) tends to a constant value 
as q — * 1~ (the calculation is detailed at the end of the 
appendix), which means 7 = 0; thus Ea. (|Bl|) holds, 
again with a = 1. The conclusion is the same for all 
B =/= J. 

Finally, we consider now the case B = J = 1. 



Then H = N(m z 



N(mi 



i 2 y )/2 ~ A/2). 



Setting e = —(1 — 8), we want to compute 
s(m x = 0,m y = 0,m z - A/2 = 1-8). Call- 
ing M = m z — A/2, we have M = (q), with 
q = cos 9 — (sin 2 9cos2<p)/2 a random variable in 
[—1, 1] with distribution p(q). The calculations in the 
next paragraph show that p(q) diverges at the boundary 
like (1 - q) 7 , with 7 = -1/4; thus Eq.fBlJ still holds, 
now with a = 3/4. 

Derivation of the exponent 7 : 

1. B = 0: we need to compute the distribution P(y = 
sin 2 9 cos 20), close to y = 1. We have, with u = cos9: 



P(y) 



1 j-2-K 

du d<p 8 D (y- (1 - w 2 ) cos 2^) , (B6) 
-1 Jo 



where 8d(x) is the Dirac delta function. Writing y — 
1—8, with 8 <C 1, we see that only the values of u such 
that u 2 < 8 contribute. Integrating over (ft we obtain: 



P(l-8) = A — 
J-^s 2 1 



du 



v - s 2((l-u 2 ) 2 -(l-8) 2 )V 2 



(B7) 



the factor of 4 in front comes from the four values of (ft 
that contribute the same amount to P. Expanding the 
denominator, neglecting order 8 2 terms and performing 
the change of variable u = s/8t, we get: 



P(l-8)^ J 



1 V2dt 



(B8) 



Since this last integral does not depend on 8 and has a 
finite value, we conclude that 7 = for B = 0. 

2. B = J = 1: we need to compute now P(y = u — 
(1/2)(1 — u 2 ) cos 2 (ft) close to y — 1. This reads: 



P(y) 



1 /.27T 

du 



y — u + -(1 — u 2 ) cos 20 



Solving for (ft inside the delta function, we get 

2(u-y) 



cos 20 



1 -u 2 



(B9) 
(B10) 
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This time, only the values of u close to u = 1 contribute 
to P; thus, we write it = 1 — s, y = 1 — 5. From the 
inequalities — 1 < cos 20 < 1, we get, neglecting terms of 
order s 2 , 8/2 < s < \/28. Integrating the delta function 
over </>, we have the expression for P: 



P(l-<5)~4 f 
Js. 



'2,S 



ds 



Is/2 ((2s-s 2 ) 2 -4(6-s) 2 ) 1/2 
the change of variable t = 2s/ 8 yields 

^ dt 



(Bll) 



P{\ - 8) 



(B12) 



This integral converges close to t = 1 ; it diverges however 
at large t, like t 1 / 2 ; since t diverges as 8~ x l 2 , we finally 
get 7 = -1/4. 



3. General case J: we do not detail here the 
calculations, which are similar to those above. As soon 
as B 7^ J, the result is 7 = 0, and thus a = 1. 
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